/*
Er moeten ook dingen wel/niet gediscount worden met rente, maar daar heb ik nog niet
echt over nagedacht. Ik weet wel zeker dat wat er nu staat niet klopt :P

Voor de rest: kunnen we de maximaal nuttige M1 en M2 niet berekenen adhv van deltaX, deltaT en S0?
De waarde van een punt verplaatst zich maximaal met 1 plaatsstap per tijdstap.
Dus later zijn punten die ver van S0 afliggen verspilde berekening.

*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <getopt.h>
#include <sys/time.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include "fd-option.h"

int main(int argc, char *argv[]) {
    fdPars p;
    fdResults r;
    struct timeval start, end;
    
    getPars(argc, argv, &p);
    
    gettimeofday(&start, NULL);
    finiteDifference(&p, &r);
    gettimeofday(&end, NULL);
    
    double d1, d2, c, delta;
    
    d1 = (log(p.S0 / p.strikePrice) + (p.r + p.sigma*p.sigma * 0.5) * p.endTime)
        / (p.sigma * sqrt(p.endTime));
    d2 = d1 - p.sigma * sqrt(p.endTime);

    if(p.digital) {
        c =  exp(-p.r * p.endTime) * gsl_cdf_ugaussian_P(d2);
        delta = exp(-p.r * p.endTime) * 1./sqrt(2. * acos(-1.)) * exp(-0.5*d2*d2) / (p.sigma * p.S0 * sqrt(p.endTime));
    } else {
        c = p.S0 * gsl_cdf_ugaussian_P(d1) - p.strikePrice * exp(-p.r * p.endTime) * gsl_cdf_ugaussian_P(d2);
        delta = gsl_cdf_ugaussian_P(d1);
    }
    double runtime = end.tv_sec - start.tv_sec +  (end.tv_usec - start.tv_usec) * 1.e-6;
    if(p.verbose) {
        printf("Time elapsed:\t\t%.3lf seconds\n", runtime);
        printf("Theoretical:\t\tvalue = %lf, delta = %lf\n", c, delta);
        printf("Simulation:\t\tvalue = %lf, delta = %lf\n", r.optionValue, r.hedge);
    } else {
        printf("%e\t%f\t%f\t%e\t%e\n", runtime, r.optionValue, r.hedge, r.optionValue - c, r.hedge - delta);
    }
    return 0;
}

void finiteDifference(fdPars *pars, fdResults *results) {
    int i, j;
    double **V;
    
    int Ntime = round( pars->endTime / pars->deltaT ) + 1;
    pars->deltaT = pars->endTime / (Ntime - 1);
    int Nspace = round( (pars->M1 + pars->M2) / pars->deltaX );        
    int xIsZeroIndex = round( pars->M1 / (pars->M1 + pars->M2) * (double)Nspace );
    
    double strikePrice = pars->strikePrice;
    
    double deltaT = pars->deltaT;
    double deltaX = pars->deltaX;
    
    double tau;
    
    double a0, a1, amin1, b0, b1, bmin1;
    
    gsl_vector *vec_a0, *vec_a1, *vec_amin1, *vec_c, *vec_nextV;
    
    double k1, k2, k3;
    double sigma = pars->sigma;
    double r = pars->r;
    
    V = (double **) malloc(Ntime * sizeof(double *));
    
    for(i=0; i < Ntime; i++) {
        V[i] = (double *) malloc(Nspace * sizeof(double));
    }
    
    // Boundary conditions for tau=0
    for(j=0; j < Nspace; j++) {
        double X = (j - xIsZeroIndex) * deltaX;
        double S = pars->S0 * exp(X);
        if(pars->digital) {
            V[0][j] = S > strikePrice ? 1 : 0;
        } else {
            V[0][j] = fmax(S - strikePrice, 0.0);
        }
    }
    
    double lowBound, upBound;
    double nextLowBound, nextUpBound;
    double maxS = pars->S0 * exp( (Nspace - xIsZeroIndex) * deltaX );
    
    k1 = (r - 0.5 * sigma * sigma) * 0.25 * deltaT / deltaX;
    k2 = 0.25 * sigma * sigma * deltaT / (deltaX * deltaX);
    k3 = -0.5 * r * deltaT;
    
    if(pars->crankNicolson == 1) { // Crank-Nicolson
        a0 = -1.0 - 2.0 * k2 + k3;
        a1 = k1+k2;
        amin1 = -k1+k2;
        
        b0 = -1.0 + 2.0 * k2 - k3;
        b1 = -k1-k2;
        bmin1 = k1-k2;
        
        
        vec_c = gsl_vector_alloc(Nspace);
        vec_nextV = gsl_vector_alloc(Nspace);
        
        vec_a0 = gsl_vector_alloc(Nspace);
        vec_a1 = gsl_vector_alloc(Nspace-1);
        vec_amin1 = gsl_vector_alloc(Nspace-1);
        
        gsl_vector_set_all(vec_a0, a0);
        gsl_vector_set_all(vec_a1, a1);
        gsl_vector_set_all(vec_amin1, amin1);
        
        for(i=0; i < Ntime-1; i++) {
            tau = i*deltaT;
            
            if(pars->digital) {
                lowBound = 0.0;
                nextLowBound = 0.0;
                upBound = exp(-r * tau);
                nextUpBound = exp(-r * deltaT) * upBound;
            } else {
                lowBound = 0.0;
                nextLowBound = 0.0;
                upBound = maxS - strikePrice * exp(-r * tau);
                nextUpBound = maxS - strikePrice * exp(-r * (tau + deltaT));
            }

            gsl_vector_set(vec_c, 0, b1*V[i][1] + b0*V[i][0] + bmin1*lowBound - amin1*nextLowBound);
            gsl_vector_set(vec_c, Nspace-1, b1*upBound + b0*V[i][Nspace-1] + bmin1*V[i][Nspace-2]- a1*nextUpBound);
            
            for(j=1; j < Nspace-1; j++) {
                gsl_vector_set(vec_c, j, b1*V[i][j+1] + b0*V[i][j] + bmin1*V[i][j-1]);
            }
            gsl_linalg_solve_tridiag(vec_a0, vec_a1, vec_amin1, vec_c, vec_nextV);
            
            for(j=0; j<Nspace; j++) {
                V[i+1][j] = gsl_vector_get(vec_nextV, j);
            }
        }
        
        gsl_vector_free(vec_a0);
        gsl_vector_free(vec_a1);
        gsl_vector_free(vec_amin1);
    } else { // FTCS
        k1 *= 2.;
        k2 *= 2.;
        k3 *= 2.;
        for(i=0; i < Ntime-1; i++) {
            tau = i * deltaT;
            
            if(pars->digital) {
                lowBound = 0.0;
                upBound = exp(-r * tau);
            } else {
                lowBound = 0.0;
                upBound = maxS - strikePrice * exp(-r * tau);
            }
            
            V[i+1][0] = (1.0 + k3) * V[i][0] + k1 * (V[i][1] - lowBound) + k2*(V[i][1] - 2.0 * V[i][0] + lowBound);
            V[i+1][Nspace-1] = (1.0 + k3) * V[i][Nspace-1] + k1*(upBound - V[i][Nspace-2]) + k2*(upBound - 2.0 * V[i][Nspace-1] + V[i][Nspace-2]);
            
            for(j=1; j < Nspace-1; j++) {
                V[i+1][j] = (1.0 + k3) * V[i][j] + k1 * (V[i][j+1] - V[i][j-1]) + k2 * (V[i][j+1] - 2.0*V[i][j] + V[i][j-1]);
            }
        }
    }

    results->optionValue = V[Ntime-1][xIsZeroIndex];
    results->hedge = 1. / (2. * pars->S0 * deltaX) * (V[Ntime-1][xIsZeroIndex+1] - V[Ntime-1][xIsZeroIndex-1]);
    
    if(pars->verbose) {
        FILE *output = fopen("results.txt", "w");
        for(j = 1; j < Nspace - 1; j++) {
            double X = (j - xIsZeroIndex) * deltaX;
            double hedge = 1. / (2. * pars->S0 * deltaX) * (V[Ntime-1][j+1] - V[Ntime-1][j-1]);
            fprintf(output, "%e\t%e\t%e\n", pars->S0 * exp(X), V[Ntime-1][j], hedge);
        }
    }
}

void getPars(int argc, char *argv[], fdPars *pars) {
    pars->crankNicolson = 0;
    pars->digital = 0;
    pars->verbose = 0;
    pars->deltaX = 1.e-2;
    pars->deltaT = 1.e-3;
    
    pars->endTime = 1.;
    pars->strikePrice = 100.;
    pars->S0 = 100.;
    pars->sigma = 0.3;
    pars->r = 0.04;
    
    pars->M1 = 1.0;
    pars->M2 = 1.0;
    
    int ch;
    while ((ch = getopt(argc, argv, "cdvn:m:t:k:S:s:r:L:U:")) != -1) {
        switch(ch) {
        case 'c': pars->crankNicolson = 1; break;   
        case 'v': pars->verbose = 1; break; 
        case 'd': pars->digital = 1; break;
        case 'n': pars->deltaX = strtod(optarg, 0); break;
        case 'm': pars->deltaT = strtod(optarg, 0); break;
        case 't': pars->endTime = strtod(optarg, 0); break;
        case 'k': pars->strikePrice = strtod(optarg, 0); break;
        case 'S': pars->S0 = strtod(optarg, 0); break;
        case 's': pars->sigma = strtod(optarg, 0); break;
        case 'r': pars->r = strtod(optarg, 0); break;
        case 'L': pars->M1 = strtod(optarg, 0); break;
        case 'U': pars->M2 = strtod(optarg, 0); break;
        }
    }
    if(pars->verbose) {
        printf("Parameters:\
                \n crankNicolson (-c)\t%d\
                \n digital (-d)    \t%d\
                \n deltaX (-n)     \t%.2e\
                \n deltaT (-m)      \t%.2e\
                \n endTime (-t)    \t%.2e\
                \n strikePrice (-k)\t%.2e\
                \n S0 (-S)         \t%.2e\
                \n sigma (-s)      \t%.2e\
                \n r (-r)          \t%.2e\
                \n M1 (-L)          \t%.2e\
                \n M2 (-U)          \t%.2e\n\n",
                pars->crankNicolson,
                pars->digital,
                pars->deltaX,
                pars->deltaT,
                pars->endTime,
                pars->strikePrice,
                pars->S0,
                pars->sigma,
                pars->r,
                pars->M1,
                pars->M2
                );
    }
}

